/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _PRESSURECORRECTION_H_
#define _PRESSURECORRECTION_H_

#include "Equations/GeneralEquation.hpp"

namespace Tuna {

template <typename T>  class PressureCorrection;
  
template <typename T>  
struct Typeinfo <PressureCorrection<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};

/*! 
 ***************************************************************************
 *  This class derives from GeneralEquation through B & N trick.
 *  In this class some particular variables and methods for scalar PDE's
 *  are defined. This equations accounts for the pressure correction used
 *  in SIMPLE-like methods for decoupling the N-S equations.
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Thu May  7 13:55:49 CDT 2009 ]
 ***************************************************************************
 */
template<typename Tscheme>
class PressureCorrection : public GeneralEquation<PressureCorrection<Tscheme > > 
{
  using GeneralEquation<PressureCorrection<Tscheme> >::name;

public:
  typedef typename Typeinfo<Tscheme>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::tiny floatTinyArray_t; 
  typedef typename TunaArray<int, Typeinfo<Tscheme>::Dim >::tiny intTinyArray_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::huge ScalarField;
  typedef SparseMatrix< Diagonal< prec_t, Typeinfo<Tscheme>::Dim > > DiagMat;

  PressureCorrection() : GeneralEquation<PressureCorrection<Tscheme> >() { }

/*!
 *  Main constructor.
 *  The scalarfield \c phi_global is a reference to a global field 
 *  variable. In this case represents the pressure 
 */
  PressureCorrection(ScalarField &phi_global,
		   DiagMat &matrix, ScalarField &b,
		   const floatTinyArray_t& deltas)
    : GeneralEquation<PressureCorrection<Tscheme> >(phi_global,matrix,b,deltas)
  {
    Dim = Typeinfo<Tscheme>::Dim;
    name = "PressureCorrection";
  }
  
  ~PressureCorrection() { }; // Nothing to do
/*!
 *  The Curiously Recursive Template Pattern (CRTP) is used.  
 *  asDerived() gets an returns a reference to the derived class.
 */
  inline Tscheme& asDerived() { return static_cast<Tscheme&>(*this); }

//  
// Set functions used to initialize the parameters.
//
  inline void setUvelocity(const ScalarField& uu){ u.reference(uu); }
  inline void setVvelocity(const ScalarField& vv){ v.reference(vv); }
  inline void setWvelocity(const ScalarField& ww){ w.reference(ww); }
  inline void setPressure(const ScalarField& pstarred){ p.reference(pstarred); }
  inline void set_du(const ScalarField& ddu) {du.reference(ddu);}
  inline void set_dv(const ScalarField& ddv) {dv.reference(ddv);}
  inline void set_dw(const ScalarField& ddw) {dw.reference(ddw);}
  
/*!
 *  Call to the functions defined in the derived classes.
 */
  inline bool calcCoefficients()
  {
    if (Dim == 2)
      return asDerived().calcCoefficients2D();
    else if (Dim == 3)
      return asDerived().calcCoefficients3D();
  }

  inline bool correction() { return asDerived().correction(); }
  inline double calcSorsum();
  inline double calcResorm();
  inline void scale(float a) { p = a * p; }
//
// Get functions
//	
  const ScalarField& getPressure() { return p; }    

  inline void printInfo() {
    std::cout << "\n +-----+"
	      << "\n | PressureCorrection< ";
    asDerived().printInfo();
    std::cout << "< " << Dim << "D  > >"
	      << "\n +-----+"
	      << "\n | Jajaja    = " 
	      << "\n +-----+";  
  } 

protected:
  int Dim;
  ScalarField u, v, w, p;
  ScalarField du, dv, dw;

};

template<typename Tscheme>
inline double PressureCorrection<Tscheme>::calcSorsum() 
{ 
  double dyz = this->dy * this->dz;
  double dxz = this->dx * this->dz; 
  double dxy = this->dx * this->dy;
  double sorsum = 0.0;
    
    if (Dim == 2)
	for(int i = this->bi; i <= this->ei; ++i ) 
	    for(int j = this->bj; j <= this->ej; ++j)
		sorsum += ( - u(i+1,j) * dyz + u(i,j) * dyz 
			    - v(i,j+1) * dxz + v(i,j) * dxz );
    
    if (Dim == 3)
	for(int i = this->bi; i <= this->ei; ++i )
	    for(int j = this->bj; j <= this->ej; ++j)
		for(int k = this->bk; k <= this->ek; ++k)
		    sorsum += ( - u(i+1,j,k) * dyz + u(i,j,k) * dyz 
				- v(i,j+1,k) * dxz + v(i,j,k) * dxz
				- w(i,j,k+1) * dxy + w(i,j,k) * dxy );
    return sorsum; 
}

template<typename Tscheme>
inline double PressureCorrection<Tscheme>::calcResorm() 
{ 
  double dyz = this->dy * this->dz;
  double dxz = this->dx * this->dz; 
  double dxy = this->dx * this->dy;
  double resorm = 0.0, smp = 0.0;
    
    if (Dim == 2)
	for(int i = this->bi; i <= this->ei; ++i ) {
	    for(int j = this->bj; j <= this->ej; ++j)
		smp = ( - u(i+1,j) * dyz + u(i,j) * dyz 
			- v(i,j+1) * dxz + v(i,j) * dxz );
	    resorm += fabs(smp);
	}
    
    
    if (Dim == 3) 
	
	for(int i = this->bi; i <= this->ei; ++i ) 
	    for(int j = this->bj; j <= this->ej; ++j)
		for(int k = this->bk; k <= this->ek; ++k) {
		    smp = ( - u(i+1,j,k) * dyz + u(i,j,k) * dyz 
			    - v(i,j+1,k) * dxz + v(i,j,k) * dxz
			    - w(i,j,k+1) * dxy + w(i,j,k) * dxy );
		    resorm += fabs(smp);
		}
    return resorm; 
}

}  // namespace Tuna


#endif //_PRESSURECORRECTION_H_













